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Abstract 

We introduce a 2-dimensional lattice model of granular matter. We use a combina- 
tion of proof and simulation to demonstrate an order /disorder phase transition in the 
model, to which we associate the granular phenomenon of random close packing. 
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0. Introduction. 

Granular materials, such as a static pile of sand or salt grains sedimented in a fluid 
such as air, exhibit interesting characteristic behavior at certain volume fractions. For 
sand in air the lowest possible volume fraction (called the random loose packing density) 
is about 0.58, and the highest possible volume fraction is about 0.74. In other words a 
sand pile can exist with volume fraction anywhere in the interval (0.58, 0.74). Within 
this range there are also: the critical state density, about 0.60, and the random close 
packing density, about 0.64 [dG]. In this paper we consider a toy model for granular 
materials, the goal being to model granular behavior near the random close packing 
density. Our results support the interpretation in [Ra] of the phenomenon of random 
close packing as an order/disorder phase transition; we show in our model that at high 
density the system is sensitive to the boundary conditions while at low density it is 
not, with a perfectly sharp transition in between. 

Our model is 2-dimensional and consists of nonoverlapping, parallel, hexagonal 
"grains" for which the centers (and corners) lie on sites of the planar triangular lattice: 
{m(l/2, y/3/2) + n(l, 0) | m, n G Z} (see Figure 1). 



Figure 1. An hexagonal grain on the triangular lattice. 

To account for the effects of gravity and friction we impose the condition that a 
conflguration is allowed or legal only if each hexagonal grain intersects one of the three 
upper edges of another hexagonal grain, such that the latter grain has a center below 
that of the former (see Figure 2). Nearest neighbor sites in the lattice have separation 
1, and the hexagons all have the same integral side length s. For the simulations 
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described below we use s = 2, 3 or 4; for our proofs any s > 1 suffices. 

We use a "grand canonical" version of the Edwards model [EO] of granular matter; 
in this version the probability Pr^ (A) of a legal configuration A of n particles in a fixed 
volume V is e^'^/Z^(/u), where G M is a parameter and Zy{ii) is the normalization 
constant (grand partition function); the "infinite volume limit" [Ru] is then taken, in 
which y ^ 

Note that this model is a variation on the hard-core lattice-gas models of classical 
statistical mechanics, introduced by Lee and Yang in [LY], which use Peierls contours 
to prove a phase transition. (See [Gi, HP] for some later developments.) Specifically, 
in this method and for "extended" hard-cores in which particles are larger than a 
single lattice site, one proves that at all sufficiently high values of fi the model exhibits 
long range positional order, being sensitive at the middle of configurations to the 
precise relative location of the distant boundary, while at all sufficiently low values of 
H the model is (easily) shown to behave as a dilute, disordered fiuid, insensitive to the 
boundary. 




Figure 2. A legal configuration (boundary hexagons are in boldface). 
For our granular model we are able to prove long range positional order for all 
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sufficiently high values of but not disordered behavior at low In place of a 
proof for the latter we have performed Monte Carlo simulations at low values of to 
demonstrate disorder. The proof for high is given in Section 1 and the numerical 
results for low are in Section 2. 

We note that there was a previous granular adaption of the old hard hexagon 
models by Monasson and Pouliquen [MP]. Their model differs in several important 
details, for instance their use of periodic boundary conditions; more important is that 
they employ their model in a study of entropy rather than random close packing. 

1. Proof for high yu. 

Consider a regular triangular lattice with distance between nearest neighbor lat- 
tice sites equal to 1. We consider configurations of hard-core parallel regular hexagons, 
where a hexagon is centered at a lattice point and has side length equal to a fixed inte- 
ger s > 1. The hexagons are all inside a square container V whose boundary consists 
of hexagons which intersect in full edges (see Figure 2). 

We call a configuration of nonoverlapping hexagons inside the boundary legal if 
each hexagon h intersects one of the three upper edges of another hexagon h' . (We 
also require that h' is centered strictly below h.) We call h' a support of h. We let 
the number of hexagons inside the boundary (called interior hexagons), n, vary, and 
fix > 0. The probability of seeing a given configuration A is Pr{A) — e^"^ jZ, where 
n is the number of interior hexagons in A and Z = Z(fi) is the normalization. (For 
simplicity the notation will ignore dependence on the container V.) 

Two hexagons h, h' are said to be linked if their intersection is a full edge (i.e. a 
line segment of length s), or if there is a sequence of hexagons ho — h,hi,..., hm — h' 
such that hi intersects /ti_|_i in a full edge for i = 0, m — 1. In particular, the 
hexagons on the boundary are all linked. We are interested in the event that the 
origin lies inside a hexagon linked to a boundary hexagon; we call this event ^lb- 

A triangle is a closed regular triangle with side length 1 and vertices at lattice 
sites. Given a configuration A inside V , we define a contour in A to be one of the 
connected components of the union of all triangles in V not covered by hexagons in 
A, and all line segments in V of length strictly less than s which are intersections of 
neighboring hexagons in A. An outer contour is a contour which intersects a boundary 
hexagon or a hexagon linked to the boundary (see Figure 3). Note that the topological 
boundary of a contour C contains a closed curve 7 which encloses an area containing 
the entire contour. We call the region enclosed by 7 the region enclosed by C. 

A sublattice is a set of points which are the centers of a collection of hexagons 
which tile the plane. There are 3s^ distinct sublattices. Note that any set of hexagons 
which are linked corresponds to a single sublattice; in particular the boundary hexagons 
define a sublattice which we call the boundary sublattice. We say that a hexagon is 
on the boundary sublattice if its center is in the sublattice defined by the boundary 
hexagons. 



3 



Definition 1. Consider a hexagon h of side lengtli s centered at tlie origin. Let S 
be a set of lattice sites in h such that 5" has exactly one representative of each 
sublattice. We are interested in the event that there is a hexagon centered in S which 
is on the boundary sublattice; we call this event Ob- 




Figure 3. An outer contour (shaded region). 

Lemma 1. If there is a hexagon centered in S which is not on the boundary sublattice, 
then this hexagon is not linked to the boundary. 

Proof. If a hexagon centered in S is linked to the boundary, then it is on the sub- 
lattice defined by the boundary, since a set of linked hexagons corresponds to a single 
sublattice. | 

Lemma 2. If there is no hexagon centered in S, or if there is a hexagon centered in 
S not linked to the boundary, then there is an outer contour C such that the origin is 
in the region enclosed by C. 

Proof. If there is no hexagon centered in S, then the origin itself is inside a contour 
and the result follows. Now assume there is a hexagon centered in S, and consider 
the contrapositive. If there is no outer contour enclosing the origin, then there is no 
contour at all enclosing the origin. Thus the hexagon centered in S is linked to the 
boundary. | 
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Corollary 1. If there is no hexagon centered in S, or if there is a hexagon centered 
in S which is not on the boundary sublattice, then there is an outer contour enclosing 
the origin. 

Proof. This follows from Lemmas 1 and 2. | 

Let [0_b]'^ be the complement of the event Ob- We will give an upper bound for 
P{[Ob]^) by using a Peierls-type argument. We will first show that the probability of 
seeing a fixed contour is exponentially small for large /U. Then we will use a counting 
argument to get an upper bound on the number of possible contours. 

The size of a contour is defined as its area in units of the area of a hexagon. We 
will see shortly that the size of a contour must be an integer. Given a contour C, 
let E be the region enclosed by C. The closures of the connected components of the 
complement of C D E will be called C -interior regions. Note that the hexagons with 
edges on the topological boundary of a C-interior region R must be linked; we say 
such hexagons are on the outside of the R, and we say the remaining hexagons in R 
are on the inside of R. 

We say that we shift a C-interior region R if we translate R while holding all other 
hexagons fixed. The translation must be given by a difference x — y where x and y are 
lattice sites. Note that the relative positions of hexagons in a C-interior region R are 
unchanged by a shift. 

Given a contour C, we say a C-interior region R is on the boundary sublattice if 
the hexagons on the outside of R are on the boundary sublattice. 

Lemma 3. For any outer contour C there is a sequence of shifts of all of the C-interior 
regions such that in the resulting configuration, each of the shifted interior regions is 
on the boundary sublattice and no hexagons overlap one another. 

Proof. For each site x on the boundary sublattice define a neighborhood of x as 
follows. Let hhe a regular hexagon of side length s centered at x. Then consists 
of all lattice sites in h except those on any of the bottom three edges of h. 

Note that the neighborhoods Nx are disjoint and together cover all the lattice sites. 
Create a sequence of shifts of the C-interior regions as follows. For each C-interior 
region i?, take a hexagon h on the outside of R; assume h is centered at y e Nx- Then 
shift Rhy X — y. Clearly, the hexagons on the outside of the shifted C-interior region 
are centered on the boundary sublattice. We must also check that the shifts do not 
create overlap. 

To this end, let hi and h2 be any two (distinct) hexagons in C-interior regions 
i?i and i?2, and let h'l and be the images of the hexagons under the shifts of Ri 
and i?2 described in the preceding paragraph. If Ri — R2 then clearly h'l and /12 do 
not overlap. Thus assume Ri ^ R2, and suppose hi and h2 are centered at yi and 
7/2, respectively. Then yi e A^xi and 1/2 £ Nx2 for some xi ^ X2, and h'l and are 
centered at Xi and X2, respectively. Since Xi and X2 are both points on the boundary 
sublattice, h'l and do not overlap, as desired. | 



5 



Note that the configuration produced by the protocol in Lemma 3 does not nec- 
essarily produce a legal configuration, just a configuration with no overlaps. 

Lemma 4. The size of a contour C is an integer. 

Proof. Consider a configuration produced by the protocol in Lemma 3. The new 
configuration has contours Ci, C2, ■ ■ ■ , Cm in place of the original contour C. Since the 
protocol creates no overlaps, and since all the shifted C-interior regions remain within 
the region enclosed by C, the contours Ci, C2, . . . , Cm have the same combined area 
as the contour C. Furthermore, the shifted C-interior regions are all on the boundary 
sublattice, so the hexagons bordering each Ci are all on the boundary sublattice. Thus 
we conclude that each Ci could be completely covered by nonoverlapping hexagons, 
all on the boundary sublattice. The result follows. | 

Lemma 5. Fix an outer contour C of size k. There is a one-to-one correspondence 
between legal configurations A with exactly n interior hexagons and C as an outer 
contour, and legal configurations A' with exactly n + k interior hexagons. 

Proof. Let Aq and Ai be two distinct configurations with the outer contour C, and 
assume Aq and Ai have uq and ni interior hexagons, respectively. Using Lemmas 3 
and 4, shift the C-interior regions of ^0 and Ai to produce configurations Aq and Ai 
which both have contours Ci, C2, . . . , Cm that can be completely covered by nonover- 
lapping hexagons. Cover these contours with nonoverlapping hexagons to produce 
configurations Aq and A'l having no + k and ni + k interior hexagons, respectively. We 
claim first that Aq and A'l are legal configurations; of course it suffices to show that 
Aq is a legal configuration. 

We have to show that each hexagon in Aq has a support. First consider a hexagon 
h in Aq in one of the shifted C-interior regions. Assume h is on the inside of the shifted 
region. Because shifts do not affect relative positions of hexagons inside the region, 
and since the configuration was legal before the shift, h must have a support. Now 
assume h is on the outside of the shifted region. If h does not have a support in the 
shifted region, then h must have had a support outside the region before shifting. Since 
the region Ci U . . . U Cm is completely filled with hexagons, one of these must be a 
support of h. Finally consider a hexagon h not in one of the shifted C-interior regions. 
Since the contours Ci, C2, . . . , Cm were completely filled with hexagons, clearly h has 
a support. 

Next, we claim that Aq and A[ are distinct configurations. Note first that the 
C-interior regions of Aq and Ai have identical outsides, because the contour C defines 
these outsides. Thus there is an obvious pairwise association between the C-interior 
regions of ^0 and the C-interior regions of Ai. Since ^0 and Ai are distinct, either 
at least one of these pairs of C-interior regions, say Rq and must have different 
insides, or ^0 and Ai must be different outside the region enclosed by C. In the latter 
case, the configurations Aq and A'l must be distinct, because the shifts done by the 
protocol in Lemma 3 do not change anything outside the region enclosed by C. In the 
former case, Rq and Ri have distinct insides. This of course does not change after 
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shifting, and so A'q and A[ are distinct. In either case Aq and A[ are distinct, so we 
have the desired correspondence. | 

Lemma 6. Let C be a fixed contour of size k. The probabihty that a configuration 
has the contour C is at most e~''^^. 

Proof. To prove this, we use the association in Lemma 5. Let Z be the normalization, 
and let Ec be the event that a configuration has the contour C. Let Hn be the number 
of legal configurations A having the contour C and n interior hexagons, and let H'^ be 
the number of legal configurations A' having n + k interior hexagons. By Lemma 5 we 
have that i/^ > -f^n, and of course we also have that X^^q ^''^^^^^H'^ < Z. Thus, we 
have the estimate 



as desired. Note that this estimate is independent of the size of the container V . | 



We are finished with half of the Peierls argument. Now we provide an upper 
bound on the number of contours of a given size. We do this by counting graphs 
whose vertices are triangles in a contour. Note that there are 6s^ triangles inside a 
hexagon. 

Before we begin the counting argument we need the following well-known facts 
from graph theory: 

Lemma 7 a. Let T be a spanning tree for a set of n points. Then T has n — 1 edges. 

b. A graph G' produced by duplicating every edge of a graph G is Eulerian. 

Now to count the contours, we make the following observation about the structure 
of a contour C. The union of all the triangles in a contour consists of several disjoint 
connected components. These components are joined to neighboring components by 
line segments in C of length strictly less than s; recall that such line segments are the 
intersections of neighboring hexagons. Thus, the minimum number of lattice segments 
in a path between triangles in neighboring components is at most s — 1, where by a 
lattice segment we mean a line segment joining nearest neighbor lattice sites. This 
leads to the following lemma. 

Lemma 8. Suppose C is a contour of size k. Then m — 6s^k is the number of triangles 
in the contour. There is a sequence (ti, . . . ,t2m-i) of triangles in C such that each 
triangle in C is some ti , and such that the mimimum number of lattice segments in a 
path joining ti and ti+i is < s — 1. 

Proof. Let Z he a set of vertices, one for each triangle in C. Define the distance 
between vertices in Z as one plus the minimum number of lattice segments in a path 
joining the corresponding triangles in C. Partition Z into Zi, . . . , Zk, where the Zj 
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correspond to the connected components of the union of aU the triangles in C. For 

each i, join two vertices in Zi by an edge iff the distance between them is 1. Then 
one by one remove edges comprising cycles in each Zj (this process is not necessarily 
unique). 

Next, for each i ^ j, join x & Zi to y & Zj by an edge iff the distance between x 
and y is less than or equal to s. By preceding considerations we see that the resulting 
graph is connected. One by one remove edges comprising cycles to produce a tree T 
spanning all the vertices of Z. Note that all (m — 1) of the edges of T have length < s. 

Now define a duplicate graph D which has the same vertices as T but which has 
two edges joining each pair of vertices which are joined by an edge in T. Then D is an 
Eulerian graph, so there is an Eulerian path, that is, a path F in D which traverses 
every edge exactly once. D has 2(m — 1) edges, so F traverses 2{m — 1) + 1 vertices, 
counting repeats. Clearly F traverses each vertex of T at least once. So take ti to be 
the triangle corresponding to the ith vertex traversed by F. | 

Lemma 9. The number of contours C of size k such that the origin is in the region 
enclosed by C is less than PmQ^, where m = 6s^k, pm = 6m^ and q = 36(s + 1)^. 

Proof. Using Lemma 8, for any contour C of size k we have a corresponding sequence 
(ti, . . .t2m-i) of triangles in C, where m = Gs'^k. Moreover, since the sequence covers 
all the triangles in C, distinct contours are associated with distinct sequences (note 
that a contour is totally defined by the positioning of its triangles). There are no more 
than 6m^ possible triangles that a contour enclosing the origin can contain, and given 
the position of the ith triangle there are at most 6(.s + 1)^ possibilities for the position 
of the {i + l)st. Since there are 2m — 1 < 2m total elements of the sequence, we may 
take q = 36(s + 1)^ and Pm — 6m^ to get the desired result. | 

Now we are ready to combine the two main ingredients of the Peierls argument 
into the final result: 

Theorem 1. The probability that there is a hexagon centered in S which is on the 
boundary sublattice goes to 1 as // goes to infinity. That is, Pr(OB) ^ 1 as /x ^ oo. 

Proof. By Corollary 1, Lemma 6 and Lemma 9, Pr{[OBY) is bounded above by 
Ylk^^i Pm{(f'^ e~^)^ , where again m = 6s^fc. Since g^* is a constant depending only on 
s, and since Pm is polynomial in k, we have that for jj. sufficiently large, the summation 
bounding Pr{\QBY) is arbitrarily small. Note that the estimate underlying this result 
is independent of the size of the container V . | 

We have abbreviated this result by saying that at sufficiently large // the system 
has long range order. We also have the following percolation result. 

Corollary 2. The probability that the origin lies inside a hexagon linked to the 
boundary goes to 1 as goes to infinity; that is, PrijdLB) — >^ 1 as /i ^ oo. 
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Proof. This follows by Lemma 2 and the same argument as in Theorem 1. | 

2. Numerical Results for low /i. 

We ran Markov chain Monte Carlo simulations on the model for a range of values 
of fi. We checked that the Monte Carlo runs were not sensitive to the initial condition 
(see Figure 4); since lower volume fraction initial conditions tended to equilibriate 
faster, we started the remaining Monte Carlo runs with void configurations. 

If /i << typical configurations do not fill the container - see Figures 5 and 6 - 
and it is harder to develop useful data. Our goal in this section is to show that in the 
infinite volume limit the boundary has no infiuence near the origin for small /x. As the 
main object of our simulation we consider the quantity defined as follows. First 
recall the set 5" defined in the preceding section, namely a set of representative lattice 
sites for each of the different sublatticcs, such that S is contained inside a regular 
hexagon of side length s centered at the origin. Recall that we define a sublattice to 
be a set of lattice sites corresponding to the centers of a collection of hexagons which 
tile the plane. 

Definition 2. For a fixed container V, we define p{iJ,) to be the probability that there 
is a hexagon h centered in S such that h is centered in the boundary sublattice. 

Note that is the same as the quantity P{Ob), but here we emphasize its depen- 
dence on fi. We want to show that in the infinite volume limit, p{p) is constant in some 
interval of positive length; its value there should be l/3s^, where is the number of 
sublattices. Recall from our results in the previous section that ]3(/x) — > 1 as /x — > cxd, 
uniformly in system size. 

Our argument will concentrate on the interval [1, 2] for Simulation for inside 
this interval and inside the interval [0, 10] suggests that in the infinite volume limit, 
p{lj,) is indeed constant inside [1,2] (see Figures 7-12). 

To obtain numerical estimates of p(//) we considered the following functions on 
our Monte Carlo runs. For a configuration A we let 5{A) = 1 if there is a hexagon 
in A centered at a point in S in the boundary sublattice; we let S{A) = otherwise. 
We define t{A) = 1 if there is a hexagon in A centered at a point in S, and t{A) = 
otherwise. 

For systems ranging in volume from 276 to 1151 (in units of hexagon volume) we 
evaluated S and t on configurations Ai, A2, . . ., and for each system size we consider 
the following statistic: 

PM := ^[S{Ai) + 5{A2) + ■■■ + d{AM)] 

where 

T = t{Ai) + t{A2) + --- + t{AM) 
So the expected value of pm is exactly p(//). We obtain confidence intervals for 
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Pm in the same way as in [AR]; in particular we determine the mixing time for our 
simulations using the biased autocorrelation function on volume fraction data, and 
then use the common method of batch means [Ge] with about 10 batches for each run, 
with batch size M chosen so that there are at least 5 mixing times per batch (except 
in the transition region). 

If the boundary has no influence near the origin, hexagons should appear in each 
sublattice with equal probability, so we expect that the limiting value of p{p) is exactly 
l/(3s^) for small /U. In Figures 7-12 we compare data from our Monte Carlo runs to 
the line y — l/(3s^), with s — 3. For fx e [0,4] the data suggests that follows 
the line; then in the range fx e [4,6], increases to about 1; for // > 6, p{fx) stays 
near the line y = 1. In Figures 11-12 we consider more detailed data for fi in [1,2]. 
Our 95% confidence intervals cover the line y = l/(3s^) more than 95% of the time, 
as appropriate. 

We note that the transition region changes as s increases. In particular, as s 
increases the smallest value of fj, such that p{iJ,) > l/(3s^) seems also to increase; 
compare Figures 9, 13, 14. 

3. Conclusion. 

Our argument is based on the behavior of p{ii) - the probability that a hexagon 
near the origin is on the same sublattice as the boundary hexagons - as a function of the 
parameter ^, the variable controlling average volume fraction. We have proven that for 
sufficiently large positive /U, p{n) is greater than l/(3s^), and in fact p{n) approaches 
1 as /X — > cxD, uniformly in the size of the system. In addition we have numerical 
evidence that in an interval above zero, p{fj,) has the constant value l/(3s^) in the 
infinite volume limit, indicative of disorder. As the two types of behavior cannot be 
connected analytically we conclude [FR] that the model undergoes a phase transition 
at some positive /i. The transition can be seen in Figure 9, but it would take much 
more simulation to demonstrate singular behavior at a specific value of Instead, our 
argument for the existence of a transition is based on failure of analyticity. (To use 
simulation to show that dependence on the boundary survives in the infinite volume 
limit requires careful study of the size of the simulation samples, while the burden is 
easier to show independence of the boundary, as we do.) 

The transition we have found is of the order/disorder type since the long range 
order which we prove to hold at large fj, is absent at low We note that, as usual 
in hard-core lattice models [HP] , our results only apply for a finite ratio s of hexagon 
size to lattice spacing; our upper bound on ordered behavior diverges as s — > oo. 
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Simulation Results 
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Figure 4. Plot of volume fraction versus number of moves, from three different initial 
volume fractions, for a system of volume 729 and ^ = 1. 
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Figure 5. A configuration of 250 hexagons in equilibrium at // = —4, in a system of 
volume 729. 
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Figure 6. A plot of a configuration in equilibrium at /U = 1 
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Figure 9. Plot of p{^) vs. n for a system of volume 276, with error bars, for s = 3. 
The hue is = ^ = j^. 
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Figure 10. Plot of p{fi) vs. /i for a system of volume 1151, with error bars, for 
The line is p{ii) = 3^ = ^■ 
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Figure 11. Plot of p{^) vs. fi for a system of volume 276, with error bars, for s = 3. 
The line is = = 
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Figure 12. Plot of vs. /i for a system of volume 1151, with error bars, for s = 3. 
The line is p{iJ,) — — 
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Figure 13. Plot of p{p) vs. n for a system of volume 276, for s = 2. The line is 
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Figure 14. Plot of vs. for a system of volume 276, for s = 4. The line is 
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